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Preservation of the maximum principle is studied for the combination of the linear finite 
element method in space and the ^-method in time for solving time dependent anisotropic 
diffusion problems. It is shown that the numerical solution satisfies a discrete maximum 
principle when all element angles of the mesh measured in the metric specified by the 
inverse of the diffusion matrix are nonobtuse and the time step size is bounded below and 
above by bounds proportional essentially to the square of the maximal element diameter. 
The lower bound requirement can be removed when a lumped mass matrix is used. In two 
dimensions, the mesh and time step conditions can be replaced by weaker Delaunay-type 
conditions. Numerical results are presented to verify the theoretical findings. 
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1 Introduction 

We are concerned with the linear finite element solution of the initial-boundary value problem (IB VP) 
of a linear diffusion equation, 

' u t ~ V • (DVu) = f(x,t), in l] T = (]x(0,T] 

< u(x,t) = g{x,t), on d(l x [0,T] (1) 

k u(x,0) = uo(x), in !] x {t = 0} 

where Q C R d (d = I, 2, or 3) is a connected polygonal or polyhedral domain, T > is a fixed time, 
f(x,t), g(x,t) and uq(x) are given functions, and O = D(a;,t) is the diffusion matrix assumed to be 
symmetric and strictly positive definite on Qt- This IB VP is a heterogeneous anisotropic diffusion 
problem when D varies from place to place and has not-all-equal eigenvalues at least on a portion of 

n T . 

Anisotropic diffusion problems arise from various areas of science and engineering including plasma 
physics [HH [201 EU EH HI] , petroleum reservoir simulation [TJ [21 El EU EH] > and image processing 
H EH2ZH320D1 EI], ibvp @ is a prototype of those anisotropic diffusion problems and satisfies the 
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maximum principle. However, when a standard numerical method such as a finite element or a finite 
difference method is used to solve this problem, the numerical solution may violate the maximum 
principe and contain spurious oscillations. It is of practical and theoretical importance to study when 
a numerical solution satisfies a discrete maximum principle (DMP) as well as develop DMP-preserving 
numerical schemes. 

The research topic has attracted considerable attention from researchers since 1970's and success 
has been made for elliptic diffusion problems; e.g, see |3lElElinil23ll25ll2Sll2BlEIll32lE3ll35l 
1361 H31 05] SHI H2 [501 [53]. For example, it is shown in [3l [8] that for isotropic diffusion problems, 
the requirement of all element angles of the mesh to be nonobtuse is sufficient for the linear finite 
element approximation to satisfy DMP. In two dimensions, this nonobtuse angle condition can be 
replaced by a weaker, so-called Delaunay condition [UJ which requires the sum of any pair of angles 
facing a common interior edge to be less than or equal to n. For anisotropic diffusion problems, 
Draganescu et al. [11] show that the nonobtuse angle condition fails to guarantee DMP satisfaction 
for a linear finite element approximation. Various techniques have been proposed to reduce spurious 
oscillations, including local matrix modification [28, 33], mesh optimization [36], and mesh adaptation 
[32 . An anisotropic nonobtuse angle condition, which uses element angles measured in the metric 
specified by D -1 instead of angles measured in the Euclidean metric (as in the nonobtuse angle 
condition), is developed in [31] to guarantee DMP satisfaction for anisotropic diffusion problems. A 
weaker, Delaunay-type mesh condition is obtained in [21] for two-dimensional problems. The results 
of [241 [3T] are extended in [35] to problems containing convection and reaction terms. 

On the other hand, less progress has been made for time-dependent problems; e.g., see [TU1 [T3| IT U 
Iia[IElini[ISl[22l[33l[3SlH2lHSlHa[521. Most of the existing research has focused on isotropic diffusion 
problems. For example, Fujii [18] considers the heat equation and shows that the time step size should 
be bounded from below and above for a linear finite element approximation to satisfy DMP when the 
mesh satisfies the nonobtuse angle condition. He also shows that the lower bound requirement can be 
removed when a lumped mass matrix is used. Thomee and Wahlbin [48j consider general anisotropic 
diffusion problems and show that a semi-discrete conventional finite element solution does not satisfy 
DMP in general. Slope limiters are employed in [12] to improve DMP satisfaction for anisotropic 
thermal conduction in magnetized plasmas. Nonlinear finite volume methods are developed by Le 
Potier [291 EH] for time dependent problems. 

The objective of this paper is to investigate conditions for the finite element approximation of IB VP 
([I]) to satisfy DMP for a general diffusion matrix function. We are particularly interested in lower 
and upper bounds on the time step size when the ^-method and the conventional linear finite element 
method are used for temporal and spatial discretization, respectively. Two types of simplicial mesh 
are considered, meshes satisfying the anisotropic nonobtuse angle condition [31J or a Delaunay-type 
mesh condition [24J. It is known that those meshes lead to DMP-satisfaction for linear finite element 
approximations to steady-state anisotropic diffusion problems. A lumped mass matrix is also studied. 
The results obtained in this paper can be viewed as a generalization of Fujii's |18j to anisotropic 
diffusion problems although such generalization is not trivial. 

The outline of this paper is as follows. In Sect. [2j the linear finite element solution of IB VP ([!]) is 
described. Sect. [3]is devoted to the development of DMP-satisfaction conditions. Numerical examples 
are presented in Sect. [4] to verify the theoretical findings. Finally, Sect. [5] contains conclusions. 
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2 Linear finite element formulation 

Consider the linear finite element solution of IB VP |T]). Assume that an affine family of simplicial 
triangulations {7~h} is given for the physical domain 0,. Define 

U g = {v€H 1 ((l) \v\ dn = g}. 

Denote the linear finite element space associated with mesh Th by Ug. A linear finite element solution 
u h (t) G Ug 1 for t G (0,T] to IB VP Q is defined by 



du h , 



n 



dx+ [ (X7v h ) T OVu h dx = [ fv h dx, \fv h G Uq 
Jci Jn 



dt 

where Uq = Ug with g = 0. This equation can be rewritten as 



2 / ( ^T vh(lx+ Yl \K\(Vv h fB K Vu h dx= f fv h dx, \/v h eU { 



(2) 



(3) 



KeT h " ^ K&T h KeT h ' 

where \K\ is the volume of element K and 

H>K = 7—r [ Odx. 
\K\ Jk 

Equation ^ can be expressed in a matrix form. Denote the numbers of the elements, vertices, and 
interior vertices of Th by N e , N v , and N v i, respectively. Assume that the vertices are ordered in such 
a way that the first N v i vertices are the interior vertices. Then Uq and u h can be expressed as 

Uq = span-Oi,--- ,<j>N vi }, 

N vi N v 

u h = J2uj(j>j+ Yl ( 4 ) 

3=1 j=N vi +l 

where cj)j is the linear basis function associated with the j' th vertex, a,j . We approximate the boundary 
condition in as 

Uj = 9j = g(dj,t), j = N vi + l,...,N v . (5) 

Substituting Q into taking v h = </>i (i = 1, ...,N v i), and combining the resulting equations with 
(ph, we obtain the linear algebraic system 

M d ^ + Au = f, (6) 
where u = {u x , ...,u Nvi ,u Nvi+1 , ...,u Nv ) T , f = (/i, f Nvi , g Nvi +i, ■■■,9n v ) T , 



M 





M 12 " 


A = 


' An 


A 12 ' 













I 



(7) 



and / is the identity matrix of size (N v — N V i). The entries of mass matrix M, stiffness matrix A, and 
right-hand-side vector / are given by 

rriij = y~*l (frjfadx, i = 1, ...,N vi , j = 1, ...,N V (8) 



KeT h 



K 



\K\ (V&) T Dif Vfa, i = l,...,N vi , j = l,...,N v (9) 



K£T h 



T. 'r-I- J K 



KeT h 
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We use the ^-method with a constant time step At for time integration. Let u n and u n+1 be the 
computed solutions at the current and next time steps, respectively. Applying the ^-method to the 
first N v i equations, we get 



[Mh M 12 



u 



n+l 



U 



At 



+ [An Ai 2 ]((l - 9)u n + 9u 



n+l^ 



where 



■n+e 



r " = [fl(tn + 9 At), f Nm (t n + 9At)] T . 

For the last N v — N v i equations (corresponding to the boundary condition) , we use 



,n+l 
l 3 



if; 1 ' = g(aj,t n+1 ), j = N vi + 1, ...,N V 



Combining (11) and (12), we have 



Bu n+1 = Cu n + At f n+d , 



where 
B = 



Mn M12 
/ 



+ 9 At 



A n A12 




c 



Mn M12 




en+e 



h(t n + 9At), ••• , f Nvi (t n + 9At), —g(a Nvi+1 ,t n+ i) 



;i -9) At 
1 



An A12 




At 



(11) 



(12) 
(13) 

(14) 
(15) 



It is worth noting that the right-hand side vector, f n+e , is formed from the values of the right-hand 
side function f(x,t) and the boundary function g(x,t). We are interested in conditions under which 
the scheme satisfies DMP. 



3 Conditions for DMP satisfaction 



In this section we develop the conditions (on the mesh and time step size) under which scheme (13) 
satisfies DMP. The main tool is a result from |46| which states that the solution of a linear algebraic 
system satisfies DMP when the corresponding coefficient matrix is an M-matrix and has nonnegative 
row sums. We first discuss the general dimensional case along with the anisotropic nonobtuse angle 
condition developed in [31] and then study the two dimensional case with the Delaunay-type mesh 
condition developed in |24j . 

We first introduce some notation. Consider a generic element K S Th and denote its vertices by 
a^, af , o,j,y Denote the face opposite to vertex af (i.e., the face not having af as its vertex) 
by Sf and its unit inward (pointing to af) normal by nf. The distance (or height) from vertex af 
to face Sf is denoted by hf. Define q- vectors as 



K 



= Tk> i = l ) ...,d+l. (16) 
hf 

Obviously, we have hf = l/\\qf\\. 

We now consider the mapping H> K ^x : K — > K; see Fig. [lj The q- vectors and heights associated 
with K are denoted by qf, and hf . We have relations 



af = ®Jaf, Sf = D~Jsf, \K\ = det(0 K )-h\K\, qf = B^qf , hf = (17) 

I 
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a 3 «i of of 

Figure 1: Sketch of coordinate transformations from K to K and to K. 



The dihedral angle between surfaces Sf and <9f (i ^ j) is denoted by a^. It can be expressed as 



cos(S#) 



I5fl 



iQi IIdkII?,- I|d k 



« / 3 



(18) 



where ||qf ||n> K = 
in the metric specified by D^ 1 . 



qf) T Ti)K(lf- Note that af- can be considered as a dihedral angle of K measured 



3.1 General dimensional case: d > 1 



We now are ready for the development of the DMP satisfaction conditions for scheme (13) for the 
general dimensional case. We first have the following four lemmas. 



Lemma 3.1 For any element K S Th and i, j = 1, d + 1, 

cos(5^ ) 



hfhf ' 
1 



for i / j 
for i = j 



(19) 



where (\>i and (f)j are the linear basis functions associated with the vertices af and , respectively. 
In two dimensions (d = 2), 

|if|(V&) T B K V<^ : 

Proof, see PI 



V / det(Dx) ..„ K . ... . . 1 _ „ 

cot(ajA i^j, i, j = 1,2,3. 



(20) 
□ 



Lemma 3.2 The stiffness matrix A defined in and |5p is an M-matrix and has nonnegative 
row sums if the mesh satisfies the anisotropic nonobtuse angle condition 



0<5£ < -, Vi,j = l,...,d + l,i^j, VKeT h . 



7T 



l 3 ~ 2 

Proof. See |3H Theorem 2.1 and its proof]. 



(21) 

□ 
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Lemma 3.3 Matrix B defined in (14) (0 < 9 < 1) is an M -matrix if the mesh satisfies (21) and 
the time step size satisfies 



At > 



1 



max max 



hfhf 



9(d + l)(d + 2) K&T h i,j=i,...,d+i cos(af-)X min (B K ) ' 



(22) 



Proof. We first show that M + 9 At A is a Z- matrix, i.e., it has positive diagonal and nonpositive 
off-diagonal entries. From Q we only need to show 



ma + 9 At an > 0, % = 1, N vi 

my + 9 At a,ij < Vi ^ j, i = 1, ...,N vi , j = 1, ...,iV„. 



(23) 
(24) 



Let Wj be the patch of the elements containing vertex aj. Notice that V<pi = when K ^ cjj. Recall 
from [7] that 

f I if I /" o . 2\K\ , . 

(25) 



K 



dx 



(d+i)(d+2)' y x 



2 if 



(d+l)(d + 2) 



Then (23) follows immediately from mti and Lemma 3.2 



For J24j), from @, ((£]), and ([25]) we have 

m ij + 9Ata ij = ^ / fafa dx + 9 At ^ |-fC| (V0j) r B# Vt/>j 



^ ( / dx + 9 At \K\ (V4>if B K V^) 
V f / da; + 9 At \K\ (V0 i/f ) T D x V^ A . 



(26) 



where ix and denote the local indices (on element K) of vertices ctj and a^. From (25) and 
Lemma 3.1 we get 



rriij + 9Ataij = \K\ 

The right-hand side term is nonpositive if 

At > 



1 - 9At C °l {a L KjKJ 

(d+l)(d + 2) hK h K 



(27) 



n , , x / , \ max max ,_ r , N 
9(d + l)(d + 2) x- 6 7h ij=i,...,<f+i cos(a£) 



hfhf 



(28) 



Moreover, (17) implies 



Thus, we have 



max 



<h¥< 



(29) 



From this, we can see that (22) implies (28). Hence, we have shown that B is a Z- matrix when (22) 
holds. 
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To show B is an M-matrix, we recall from ( 14 ) that 



B 



Mn + 6 At An Mii + 9AtAi2 
/ 



The fact that B is a Z-matrix means that Mn + 9 At An is also a Z-matrix and Myi + 9AtAyi < 0. 
It is easy to show that Mn + 9 At An is positive definite, which in turn implies Mn + 9 At An is an 
M-matrix. Notice 



B- l = 

This means B^ 1 > and hence B is an M-matrix 



(Mn + 9AtAii)~ 1 -{Mn + 9AtAn)- 1 {Mi2 + 9AtA u ) 
I 



Lemma 3.4 Matrix C defined in (14) (0 < 9 < 1) is nonnegative if the mesh satisfies (21) and 
the time step size satisfies 



At < 



mm mm 



(h 



K\2 



(1 -6){d+ l)(d + 2) KeT h i=l,...,d+l Xmax(P K ) ' 



(30) 



Proof. For off-diagonal entries (i 7^ j, i = l,...,N v i, j = 1,...,N V ), mij — (1 — 9)Ataij, are 
nonnegative since a^- < under condition (21) (cf. Lemma 3.2) and mij > from definition Q. To 
see if the diagonal entries are also nonnegative, from Q, (|9), and (25) we have 

2 (l-0)At\ 



Tfli 



(1 - 9) At at 



Em 



(d+l)(d + 2) ( h K)2 J' 



(31) 



The right-hand side term is nonnegative if 



At < 



min min (h 



K\2 



;i-0)(d+l)(d + 2) xer ft i=i,...,d+i 



From (29) we see that this condition holds when (30) is satisfied. 



We are now in a position to prove our first main theoretical result. 



Theorem 3.1 Scheme ( 13) satisfies a discrete maximum principle if the mesh satisfies the anisotropic 



nonobtuse angle condition (21) and the time step size satisfies (22) and (30), i.e., 

At 



max max 



hfhf 



9(d + l)(d + 2) KeT h i.j [....,!■ \ cos(5^)A. 



< 



¥=3 

2 



minWK 



<■ mm mm 

- (1 - 9){d+ l)(d + 2) Ker ft i=i,..,d+i A ma;c (D x ) 



Proof. Scheme ( 13 ) can be expressed as 



(32) 



AU = F, 



(33) 
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where 





"JO 




r i/° " 




Uq 






-C B 




u 1 




Atf 




A = 


-C B 


U = 


u 2 


F = 


Atf 1+e 


(34) 




-C B _ 








Atf N-i+e 





and B and C are defined in (14). Scheme (33) satisfies a DMP if coefficient matrix A is an M-matrix 



and has nonnegative row sums. From Lemmas 3.3 and 3.4 we know that B is an M-matrix and C > 0. 
As a result, A is a Z-matrix. Moreover, we can show A -1 > 0. Indeed, from (33) we know that uq > 
implies u° > 0. Next, from the scheme we have u 1 = B^Atf + B^Cu . Recall that C > and 
B is an M-matrix and thus -B -1 > 0. Combining these results, we can conclude that f > implies 
u > 0. Similarly, we can show u n > if f l - l+e > 0, n = 2, ...,N. Thus, we have shown that F > 
implies U > 0. This implies A -1 > and A is an M-matrix. 

We notice that the sum of each of the second to the last (block) rows is 



B-C 



At A n 




At A 12 
I 



Since A has nonnegative row sums (cf. Lemma 3.2), A has nonnegative row sums. As a consequence, 
scheme (33) (and hence (13)) satisfies DMP. □ 



Remark 3.1 It is known (31] that a mesh, generated as a uniform mesh in the metric specified 
by M K = Ok^k 1 for all K eT h , where 9 K is an arbitrary piecewise constant, scalar function defined 
on O, satisfies the anisotropic nonobtuse angle condition (21). The reader is referred to |31j for more 
information on the generation of such meshes. □ 



The lower bound requirement on At in (32) can be avoided by using a lumped mass matrix. In this 



case, scheme (13) is modified into 



Mn 
/ 



+ 6 At 



Au A u 




u 



n+l 



Mn 




An A X2 




(\-e)At 

where Mn is the lumped mass matrix with diagonal entries 



u n + Atf n+e , 



(35) 



l,...,N vi . 



The following theorem can be proven in a similar manner as for Theorem 3.1 



Theorem 3.2 Scheme (35) with a lumped mass matrix satisfies a discrete maximum principle if 

(36) 



the mesh satisfies the anisotropic nonobtuse angle condition \21y and the time step size satisfies 

1 • ■ {hf? 



At < 



mm mm 



(1 - 0)(d + 1) K&T h i=l,...,d+l Xmax(PK) 
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3.2 Two dimensional case: d = 2 

We now consider the two dimensional case. It is known [23] that a Delaunay-type mesh condition, 
which is weaker than the nonobtuse angle condition (21), is sufficient for a linear finite element 
approximation to satisfy DMP in two dimensions for steady-state problems. It is interesting to know 
if this is also true for time-dependent problems. 

Consider an arbitary interior edge e^. Denote the two vertices of the edge by and aj and the 
two elements sharing this common edge by K and K' . Let the local indices of the vertices on K be ix 
and jk- The angle of K opposite eij is denoted by af K * (when measured in the Euclidean metric) 
and by „• when measured in the metric O"^ 1 . Similarly, we have a?' „• and a? „• . 

■> l K,]K A l K''JK' l K'iJK' 

Lemma 3.5 The stiffness matrix A defined in £?p and |5p is an M-matrix and has nonnegative 
row sums if the mesh satisfies the Delaunay-type mesh condition 



a ? K ,3K + arccot 



~K' 

+ a,- , , + arccot 



det(Djf) 



7) COt Kx,jJ 




det(D^) 



cot(af , ,• " 



< 7r, V interior edges eu. (37) 



Proof. See [231 Theorem 4.1]. 



□ 



Lemma 3.6 Matrix B defined in (14) (0 < 9 < 1) is an M-matrix if the mesh satisfies (31) and 
the time step size satisfies 

2 \K\ + \K'\ 



At > 



max ■ 



(38) 



0(d+ l)(d + 2) ey v /det(D A ')cot(5^ jK ) + s/det(p K >) cot(ag #Jjc/ ) ' 

where the maximum is taken over all interior edges and K and K' are the two elements sharing the 
common edge e^j. 



Proof. Inequality (38) follows from (20), d25b, and &261. 



□ 



Lemma 3.7 Matrix C defined in (14) (0 < < 1) is nonnegative if the mesh satisfies (31) and 
the time step size satisfies 

2 



At < 



;i-0)(d+l)(d + 2) 



mm 



E |^|A ma ,(Dx 



(39) 



where the minimum is taken over all interior vertices and uji is the patch of the elements containing 
dj as its vertex. 



Proof. The proof is similar to that of Lemma 3.4 Indeed, Lemma 3.5 implies that the off-diagonal 
entries of C are nonnegative under condition (37). For diagonal entries, from (31) we get 

I if I 



(1 - 0)At ai 



2\uji 



(d+l)(d + 2) 



{I -9) At 



From (29), we can see that the right-side term of the above equation is nonnegative when (39) holds. 

□ 



Using the above results we can prove the following theorems in a similar manner as for Theorems 3.1 
andE21 
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Theorem 3.3 In two dimensions, scheme (13) satisfies a discrete maximum principle if the mesh 



satisfies the Delaunay-type mesh condition (31) and the time step size satisfies (38) and (39), i.e., 
2 \K\ + \K'\ 



max ■ 



6{d + l)(d + 2) ey v / det(D x )cot(5^ ii/f ) + v / det(D x cot (5^ 

2 



< At < 



mm 



(l-0)(d+l)(d + 2) i £ |-K"| WOOk) (/»: 



r, (40) 



where the maximum is taken over all interior edges, K and K' are the two elements sharing the 
common edge eij, and the minimum is taken over all interior vertices and uji is the patch of the 
elements containing a,i as its vertex. 



Theorem 3.4 In two dimensions, scheme (35) with a lumped mass matrix satisfies a discrete 



maximum principle if the mesh satisfies the Delaunay-type mesh condition (31) and the time step size 
satisfies 



At < 



[i-e)(d + i; 



mm 



£ \K\\ max (B K )(h? K )-2' 



(41) 



4 Numerical results 

In this section we present numerical results obtained for two examples in two dimensions to demon- 



strate the significance of both mesh conditions (21) and (37) and time step conditions (40) and (41) 



for DMP satisfaction. Three types of mesh are considered. The first is denoted by Mesh45 where 
the elements are isosceles right triangles with longest sides in the northeast direction. The second 
one is denoted by Meshl35 where the elements are isosceles right triangles with longest sides in the 
northwest direction. Examples of Mesh45 and Meshl35 are shown in Figs, ga) and (b). The third 
type of mesh, denoted by MdmPi is a uniform mesh in the metric Mfjmp = ^ which guarantees 



satisfaction of mesh condition (|21[) (cf. Remark 3.1). 



The implicit Euler method (corresponding to 9 = 1 in (13)) is used in our computation. For this 



method, conditions (32), {[361) , (40), and (41) place no constraint on the upper bound of At. For this 



reason, we consider only the lower bound for the time step size. The lower bound in (32) (related 



to the anisotropic nonobtuse angle condition) is denoted by AtAni and that in (40) (related to the 
Delaunay-type mesh condition) by At Del- Unless stated otherwise, the presented results are obtained 
after 10 steps of time integration. 

Example 4.1 The first example is in the form of IB VP with 

f = 0, fi = [0,1] 2 \[0.4,0.6] 2 , 9 = 0onr mt , g = 4 on T in , 

where T ou t and Tj n are the outer and inner boundaries of O, respectively; see Fig. |3^a). The initial 
solution Uq(x, y) is given as 



uo(x,y) 



4, on T in 

0, in O\[0.2,0.8] 2 

increases linearly, from T m jd to Tj 
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(a): Mesh45 (b): Meshl35 

Figure 2: Examples of Mesh45 and Meshl35. 



where r mi ^ is the boundary of subdomain [0.2, 0.8] 2 ; see Fig. [3J The diffusion matrix is taken as 



50.5 49.5 
49.5 50.5 

which has eigenvalues 100 and 1. The principal eigenvectors are in the northeast direction. 

This example satisfies the maximum principle and the exact solution (whose analytical expression 
is unavailable) stays between and 4. Our goal is to produce a numerical solution which also satisfies 
DMP and stays between and 4. 

We first consider Mesh45 and Meshl35. Mesh45 satisfies the anisotropic nonobtuse angle condition 



(21) since its maximum angle in the metric M 



-1 : 



is 0.477T. It is known [24J that (21) implies the 



Delaunay-type mesh condition, (37). Indeed, by direct calculation we can find that the maximum of 



the left-hand-side term of (37) is 0.947T. On the other hand, Meshl35 satisfies neither of (21) and (37), 



with the maximum angle in the metric M = D 1 being 0.94-7T and the maximum of the left-hand-side 



term of (37) being 1.877T. 



The solution contours (after 10 time steps) using Mesh45 and Meshl35 with h = 2.5 x 10 -2 and 
At = 1.5 x 10 -4 are shown in Fig. [4J where h denotes the maximal height of triangular elements of 
the mesh and u m i n is the minimum of the numerical solution. No undershoot occurs in the numerical 
solution obtained with Mesh45. 

The results for Mesh45 are listed in Table [lj They show that for meshes with h < 2.5 x 10 -2 , 
At D e i is smaller than the step size At = 1.5 x 10 -4 used in the computation. As a consequence, 
time condition (40) (and mesh condition (|37[)) is satisfied and Theorem 3.3 implies that the numerical 
solution satisfies DMP. Table [T] confirms that no undershoot occurs in the numerical solution or 
Umin = 0. On the other hand, for h = 5.0 x 10 -2 , neither of time conditions (32) and (40) is satisfied 
and undershoot with u m i n = — 1.41 x 10~ 7 is observed. 

The table also records the numerical results obtained for h = 2.5 x 10 -2 and h = 1.25 x 10 -2 with 
decreasing At. One can see that no undershoot occurs when At > Atjj e i. However, undershoot occurs 
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u = 

(a): Physical domain and boundary (b): Initial solution 

Figure 3: The physical domain, boundary condition, and initial solution for Example 4.1 



when At continues to decrease and pass AtDel- This is consistent with Theorem 3.3 

It is pointed out that At Del < AtAni for all the cases listed in the table. Moreover, for some cases 
we have At Del < At < AtAni and no undershoot occurs in the numerical solution. These indicate 



that time condition (40) (related to the Delaunay-type mesh condition) is weaker than (32) (related 



to the anisotropic nonobtuse angle condition) 



Recall that Meshl35 does not satisfy mesh condition (21 ) nor (37). Thus, there is no guarantee that 



the numerical solution obtained with Meshl35 satisfies DMP. Indeed, Table [2] shows that undershoot 
occurs in all numerical solutions obtained with various sizes of Meshl35 and various At. 

Next we consider Mdmp meshes which are generated as (quasi-) uniform ones in the metric specified 
by M = ID 



Recall from Remark 



3.1 



that such meshes satisfy the anisotropic nonobtuse angle condi- 
tion (21). In our computation, Mdmp meshes are generated using BAMG (bidimensional anisotropic 
mesh generator) code developed by Hecht |23j . An example is shown in Fig. [5^b). Notice that the 
elements are aligned with the principal diffusion direction (northeast). Since the diffusion tensor D 
is constant, the mesh is generated initially based on Mdmp = ^>~ 1 and then kept for the subsequent 
time steps. 

The results obtained with Mdmp meshes are similar to those obtained with Mesh45. For example, 
for the Mdmp mesh shown in Fig. [5] (b), it is found numerically that AtAni = 4.30 x 10" 2 and 



At D el = 1-63 x 10~ d . Theorem 3.3 ensures that no undershoot occurs in the numerical solution when 
At > AtDel- R is emphasized that ( |37[ ) and (40) are not necessary for DMP satisfaction and the 



numerical solution may be free of undershoot for some smaller values of At. In fact, no undershoot 
is observed numerically for At > 10~ 4 . A undershoot-free solution obtained with the mesh shown 
in Fig. [5] (b) and time step size At = 1.5 x 10 -4 is shown in Fig. [5] (a). For the same mesh with 
At = 1.0 x 10 ^, undershoot is observed with u m i n = -1.45 x 10" b . 

Finally, we consider the lumped mass method. Theorem [3]4] implies that there is no constraint placed 
on At for the DMP satisfaction of the numerical solution with the lumped mass matrix and implicit 
Euler discretization. Indeed, for all Mesh45 meshes and At considered in Table [TJ no undershoot is 
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(a): Mesh45, u min = (b): Meshl35, u m %n = -6.57 x 1(T 2 

Figure 4: Solution contours obtained for Mesh45 and Meshl35 with h = 2.5 x 10 -2 and At = 1.5 x 10 



for Example 



4.1 



observed numerically for the lumped mass method. The same also holds for Mdmp meshes. For 
example, for the mesh shown in Fig. [5^b), no undershoot is observed in the numerical solution for 
At = 10 -4 , 10 -5 , and 10 -6 . For Meshl35 meshes, mesh condition (21) or (37) is not satisfied. 
Thus, Theorem 3.4 does not hold. For example, for a Meshl35 mesh with h = 1.25 x 10 -2 and 
the numerical solution violates DMP and has a minimum u m i n = —1.60 x 10 -2 . 



At 



Theorem 
1.5 x 10~ 4 



Example 4.2 The second example is the same as Example 4.1 except that the diffusion matrix 
is taken as a function of x and y, i.e., 



cos ( 
sin ( 



— sm ( 
cos 9 



100 
1 



cos t> 
— sin( 



sinf 
cos I 



where 6 = 6(x, y) is the angle of the tangential direction at point (x, y) along circles centered at 
(0.5, 0.5). This diffusion matrix D also has eigenvalues 1 and 100 but has its principal eigen-direction 
along the tangential direction of circles centered at (0.5, 0.5). A physical example with such a diffusion 
matrix is the toroidal magnetic field in a Tokamak device confining fusion plasma |41j . This problem 
also satisfies the maximum principle and the solution stays between and 4. 

For this example, neither Mesh45 nor Meshl35 (cf. Fig. [2]) satisfies the anisotropic nonobtuse angle 
condition (37). In the metric specified by M = D _1 , the maximum of the left-hand side of inequality 
(37) is 1.87-7T for both Mesh45 and Meshl35. Due to the symmetry of the diffusion matrix, both 



Mesh45 and Meshl35 lead to almost the same results for this example except that undershoot occurs 
at different locations. Fig. [6] shows the results obtained with these meshes for At = 5 x 10~ 5 . 

Table [3] list numerical results obtained with Mesh45 and Mdmp meshes. Recall that Theorem 3J3 
does not apply to Mesh45 meshes since they do not satisfy (37). As a matter of fact, numerical 



solutions obtained with this type of meshes with or without mass lumping violate DMP and exhibit 
undershoot. On the other hand, Mdmp meshes generated with M = D _1 satisfy the mesh condition. 
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Table 1: Numerical results obtained with Mesh45 for Example 



4.1 



h 


At^ni 


At_D e Z 


At 




5.0e-2 


1.48e-3 


3.79e-4 


1.5e-4 


-1.41e-7 


2.5e-2 


3.70e-4 


9.47e-5 


1.5e-4 





1.25e-2 


9.25e-5 


2.37e-5 


1.5e-4 





6 25e-3 


2.31e-5 


5 92e-6 


1.5e-4 


o 


3.125e-3 


5.78e-6 


1.48e-6 


1.5e-4 





2.5e-2 


3.70e-4 


9.47e-5 


1.5e-4 





2.5e-2 


3.70e-4 


9.47e-5 


1.0e-4 





2.5e-2 


3.70e-4 


9.47e-5 


5.0e-5 


-7.91e-10 


1.25e-2 


9.25e-5 


2.37e-5 


1.5e-4 





1.25e-2 


9.25e-5 


2.37e-5 


1.0e-5 


-1.31e-6 



Table 2: Numerical results obtained with Meshl35 for Example 



4.1 



h 


At^ni 




At 




5.0e-2 


1.48e-4 


2.08e-6 


1.5e-4 


-8.99e-2 


2.5e-2 


3.70e-5 


5.21e-7 


1.5e-4 


-6.57e-2 


1.25e-2 


9.25e-6 


1.30e-7 


1.5e-4 


-1.58e-2 


1.25e-2 


9.25e-6 


1.30e-7 


1.0e-7 


-2.26e-2 


6.25e-3 


2.31e-6 


3.26e-8 


5.0e-4 


-1.59e-3 


6.25e-3 


2.31e-6 


3.26e-8 


1.5e-5 


-1.43e-2 


6.25e-3 


2.31e-6 


3.26e-8 


1.5e-6 


-2.11e-2 



For the lumped mass method, no undershoot occurs in the numerical solution for all values of At. 



This is consistent with Theorem 3.4. For the standard finite element method, there is no undershoot 
for relatively large At. It is interesting to point out that for this example with variable D, the lower 
bounds At^ni and Atp^ are far too pessimistic. A several magnitude smaller At can still lead to 
numerical solutions free of undershoot. 



5 Conclusions 



In the previous sections we have studied the conditions under which a full discretization for IB VP ([I]) 
with a general diffusion matrix function satisfies a discrete maximum principle. The discretization 
is realized using the ^-method in time and the linear finite element method in space. The main 



theoretical results are given in Theorems 3.1 3.2, 3.3 and 3.4 



Specifically, the numerical solution obtained with the full discrete scheme satisfies a discrete max- 



imum principle when the mesh satisfies the anisotropic nonobtuse angle condition (21) and the time 



step size satisfies condition (32). As shown in [31], a mesh satisfying (21) can be generated as a 
uniform mesh in the metric specified by a D _1 with a being a scalar function defined on FIt- On the 



other hand, condition (32) essentially requires the time step size to satisfy 



Cih 2 < At < 



Co 



(42) 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(a): u r 







Figure 5: An M DM p mesh (with N e = 2362 and N v 



with At = 1.5 x 10 4 for Example 



4.1 



(b): M DM p mesh, N e = 2362 
1357) and the corresponding solution obtained 



where C\ and C2 are positive constants, h is the maximal element diameter, and 9 £ (0, 1] is the 
parameter used in the ^-method. Obviously, this condition is restrictive. This is especially true when 
the numerical scheme with 6 £ [0.5,1] is known to be unconditionally stable and no constraint is 



placed on At for the sake of stability. Moreover, (32) indicates that the numerical solution will violate 
the maximum principle as At — > 0. This is consistent with the finding of Thomee and Wahlbin [48J 
that a semi-discrete standard Galerkin finite element solution violates DMP since the semi-discrete 
scheme can be considered as the limit of the full discrete scheme as At —> 0. Furthermore, Theorems 



3.2 and 3.4 show that the lower bound requirement on At can be removed when a lumped mass matrix 



is used. Finally, in two dimensions, the mesh and time step conditions can be replaced with weaker 
conditions (37) and (40), respectively. Numerical results in Sect. E] confirm the theoretical findings. 
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